library(tidyverse)
library(openxlsx) #read excel file
library(ggrepel) #label text
library(vegan) #Biodiversity Index
library(sf) #Map
library(scales)
library(treemapify)
source("Custom_function.R")
Use read.xlsx() function to load data from Excel file to
R.
#library(openxlsx)
filename <- "Plot_NC-BD_20230925_v4.xlsx"
# Bidoup plot
BDplot <- read.xlsx(filename, sheet = "BD_census2")
This dataset is 2 hectares:
BD-1 (1ha) is concentrated (100x100m)
BD-2 (1ha) is taken from 25 scattered subplots (20x20m).
Use n_distinct() to counts the number of distinct value
in column: family, genus, scientificName. And
group_by() to count at each plot.
BDplot |>
group_by(PlotID) |>
summarise("family" = n_distinct(family),
"genus" = n_distinct(genus),
"species" = n_distinct(scientificName),
"stem" = n())
## # A tibble: 2 × 5
## PlotID family genus species stem
## <chr> <int> <int> <int> <int>
## 1 BD-1 46 70 124 12207
## 2 BD-2 53 94 169 11575
We usually group trees into 3 size-class of DBH:
small tree [1-5)cm;
medium tree [5-20)cm;
big tree >=20 cm
BDplot <- BDplot |>
mutate(
#convert dbh from mm to cm
dbh_cm = dbh / 10,
gDBH = case_when(dbh_cm >= 20 ~ ">=20",
dbh_cm >= 5 ~ "[5-20)",
TRUE ~ "[1-5)")
)
head(BDplot)
We count the number of stem at each diameter class and each plot
(stemByDBH <- BDplot |>
group_by(PlotID, gDBH) |>
summarise(n = n()) |>
mutate(freq = round(n / sum(n), 3)) |>
arrange(desc(n)))
## # A tibble: 6 × 4
## # Groups: PlotID [2]
## PlotID gDBH n freq
## <chr> <chr> <int> <dbl>
## 1 BD-1 [1-5) 9736 0.798
## 2 BD-2 [1-5) 9478 0.819
## 3 BD-1 [5-20) 2057 0.169
## 4 BD-2 [5-20) 1749 0.151
## 5 BD-1 >=20 414 0.034
## 6 BD-2 >=20 348 0.03
Figure
ggplot(BDplot, aes(factor(gDBH), fill = factor(PlotID))) +
geom_bar(position = "dodge") +
theme_bw(14) +
labs(x = "DBH (cm)", y = "Nb of stem", fill = "PlotID")
We can also count the number of species at each diameter class.
(
SpByDBH <- BDplot |>
group_by(PlotID, gDBH) |>
summarise(
"family" = n_distinct(family),
"genus" = n_distinct(genus),
"species" = n_distinct(scientificName)
)
)
## # A tibble: 6 × 5
## # Groups: PlotID [2]
## PlotID gDBH family genus species
## <chr> <chr> <int> <int> <int>
## 1 BD-1 >=20 28 35 53
## 2 BD-1 [1-5) 45 69 120
## 3 BD-1 [5-20) 42 61 96
## 4 BD-2 >=20 33 45 64
## 5 BD-2 [1-5) 53 91 160
## 6 BD-2 [5-20) 48 74 121
Figure
Comparison of family numbers between two plots
SpByDBH |>
ggplot(aes(x = gDBH, y = family, fill = PlotID)) +
geom_col(position = position_dodge()) +
geom_text(aes(x = gDBH, y = family, label = family),
position = position_dodge(width = 0.9)) +
theme_bw(15) +
labs(x = "DBH class", y = "Number", fill = "Plot")
Comparison of genus numbers between two plots
SpByDBH |>
ggplot(aes(x = gDBH, y = genus, fill = PlotID)) +
geom_col(position = position_dodge()) +
geom_text(aes(x = gDBH, y = genus, label = genus),
position = position_dodge(width = 0.9)) +
theme_bw(15) +
labs(x = "DBH class", y = "Number", fill = "Plot")
Comparison of species numbers between two plots
SpByDBH |>
ggplot(aes(x = gDBH, y = species, fill = PlotID)) +
geom_col(position = position_dodge()) +
geom_text(aes(x = gDBH, y = species, label = species),
position = position_dodge(width = 0.9)) +
theme_bw(15) +
labs(x = "DBH class", y = "Number", fill = "Plot")
BDplot.BD1 <- BDplot |> filter(PlotID == "BD-1")
Draw a blank plot with draw_blank_plot_20ha()
function.
p <- draw_blank_plot_20ha()
p
Add add tree to map
p + geom_point(data = BDplot.BD1,
aes(x = gx, y = gy, size = dbh_cm, color = gDBH),
alpha = 0.7) +
xlim(80, 180) + ylim(400, 500) +
labs(color = "DBH class", size = "DBH(cm)")
For big tree
BDplotA <- BDplot.BD1 |> filter(gDBH == ">=20")
p + geom_point(data = BDplotA,
aes(x = gx, y = gy, size = dbh_cm),
alpha = 0.7, shape = 21,fill = "blue", color = "white") +
xlim(80, 180) + ylim(400, 500) +
labs(size = "DBH(cm)")
For medium tree
BDplotB <- BDplot.BD1 |> filter(gDBH == "[5-20)")
p + geom_point(data = BDplotB,
aes(x = gx, y = gy, size = dbh_cm),
alpha = 0.7, shape = 21,fill = "blue", color = "white") +
xlim(80, 180) + ylim(400, 500) +
labs(size = "DBH(cm)")
For small tree
BDplotC <- BDplot.BD1 |> filter(gDBH == "[1-5)")
p + geom_point(data = BDplotC,
aes(x = gx, y = gy, size = dbh_cm),
alpha = 0.7, shape = 21,fill = "blue", color = "white") +
xlim(80, 180) + ylim(400, 500) +
labs(size = "DBH(cm)")
For a particular species, px: Pinus krempfii
## So do phan bo toan bo cay trong o mau
BDplotSP <- BDplot.BD1 |> filter(scientificName == "Pinus krempfii")
p + geom_point(data = BDplotSP,
aes(x = gx, y = gy, size = dbh_cm, color = gDBH),
alpha = 0.7) +
xlim(80, 180) + ylim(400, 500) +
labs(color = "DBH class", size = "DBH(cm)")
We can also try with a detail DBH class
BDplot.BD1 <- BDplot.BD1 |>
mutate(
dbh_cm = dbh / 10 , # dbh in mm
gDBH2 = case_when(
dbh_cm >= 100 ~ ">=100",
dbh_cm >= 90 ~ "[90-100)",
dbh_cm >= 80 ~ "[80-90)",
dbh_cm >= 70 ~ "[70-80)",
dbh_cm >= 60 ~ "[60-70)",
dbh_cm >= 50 ~ "[50-60)",
dbh_cm >= 40 ~ "[40-50)",
dbh_cm >= 30 ~ "[30-40)",
dbh_cm >= 20 ~ "[20-30)",
dbh_cm >= 10 ~ "[10-20)",
TRUE ~ "[1-10)")
)
# Count number of stem in each DBH class
stemByDBH2 <- BDplot.BD1 |>
count(gDBH2) |>
mutate(freq = round(n / sum(n), 3)) |>
arrange(desc(n))
stemByDBH2
## gDBH2 n freq
## 1 [1-10) 11169 0.915
## 2 [10-20) 624 0.051
## 3 [20-30) 246 0.020
## 4 [30-40) 123 0.010
## 5 [40-50) 35 0.003
## 6 [70-80) 3 0.000
## 7 >=100 2 0.000
## 8 [50-60) 2 0.000
## 9 [60-70) 1 0.000
## 10 [80-90) 1 0.000
## 11 [90-100) 1 0.000
# Figure
stemByDBH2 |>
ggplot(aes(x = gDBH2, y = n)) +
geom_col() +
geom_text(aes(label = n), vjust = -0.5) +
theme_bw(14) +
labs(x = "DBHClass (cm)", y = "Number of stem")
The Importance Value Index (IVI) in Ecology is the quantitative measure of how dominant a species is in a given ecosystem. IVI is calculated from 3 values:
relative density of each species (RD)
relative basal area of each species (RBA)
relative frequency of each species (this can be omitted in 1 plot)
\[IVI={\left(relative.density + relative.basal.area \right)}\]
Basal area (BA) of each species can be calculated from diameter measured: \(BA = pi * R^2\)
Relative basal area (RBA): \(RBA = BA / sum(BA)\)
Relative density (RD) = number of individuals of each species / total
BDplot <- BDplot |>
mutate(BA = ((dbh/2)/1000)^2*pi) # unit: m2)
BDplot.BD1 <- BDplot |> filter(PlotID == "BD-1")
BDplot.BD2 <- BDplot |> filter(PlotID == "BD-2")
IVI.1 <- calcIVI(BDplot.BD1)
head(IVI.1)
## # A tibble: 6 × 7
## scientificName n sG UT MD IVI rank
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Lithocarpus pierrei 966 6.74 12.7 7.91 20.6 1
## 2 Polyosma nhatrangensis 1175 2.33 4.40 9.63 14.0 2
## 3 Syzygium rubicundum 726 3.21 6.06 5.95 12.0 3
## 4 Pinus krempfii 89 5.28 9.96 0.729 10.7 4
## 5 Castanopsis echinocarpa 539 3.05 5.75 4.42 10.2 5
## 6 Lasianthus saprosmoides 1076 0.236 0.444 8.81 9.26 6
IVI.2 <- calcIVI(BDplot.BD2)
head(IVI.2)
## # A tibble: 6 × 7
## scientificName n sG UT MD IVI rank
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Lithocarpus coinhensis 357 3.94 7.84 3.08 10.9 1
## 2 Castanopsis echinocarpa 483 3.27 6.51 4.17 10.7 2
## 3 Polyosma nhatrangensis 705 1.31 2.60 6.09 8.69 3
## 4 notdet 883 0.127 0.253 7.63 7.88 4
## 5 Rhaphiolepis poilanei 436 1.07 2.13 3.77 5.90 5
## 6 Lithocarpus pierrei 262 1.40 2.79 2.26 5.05 6
Figure
figureIVI(IVI.1)
figureIVI(IVI.2)
biodivIndex(IVI.1)
## Richness Simpson Shannon Evenness
## 1 124 0.9598452 3.765398 0.7811574
biodivIndex(IVI.2)
## Richness Simpson Shannon Evenness
## 1 169 0.9756952 4.200136 0.8187561
Above-ground biomass (AGB)
Below-ground biomass (BGB)
Chave et al (2005). Tree allometry and improved estimation of carbon stocks and balance in tropical forests. Oecologia, 145(1): 87–99.
AGB=Total aboveground biomass (kg)
D=DBH (cm)
H=Height (m)
\(ρ\)=Wood density (g/cm3)
BA=Basal area (cm2)
| Habitat | Equation in Chave et al (2005) |
| Tropical forest DRY with Height | \(AGB = exp(-2.187 + 0.916 * ln(ρ * D^2 * H)\) |
| Tropical forest DRY only DBH | \(AGB = ρ * exp(-0.667 + 1.784 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)\) |
| Tropical forest MOIST with Height | \(AGB = exp(-2.977 + ln(ρ * D^2 * H))\) |
| Tropical forest MOIST only DBH | \(AGB = ρ * exp(-1.499 + 2.148 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)\) |
library(BIOMASS)
woodensity <- getWoodDensity(
genus = BDplot.BD1$genus,
species = BDplot.BD1$scientificName,
family = BDplot.BD1$family
)
head(woodensity)
BDplot.BD1$meanWD <- woodensity$meanWD
Bidoup Plot: Tropical forest MOIST only DBH
\(AGB = ρ * exp(-1.499 + 2.148 * ln(D) + 0.207 * (ln(D))^2 – 0.0281 * (ln(D))^3)\)
BDplot.BD1 <- BDplot.BD1 %>%
mutate(dbh_cm = dbh/10,
AGB = meanWD * exp(-1.499 + 2.148*log(dbh_cm) + 0.207*log(dbh_cm)^2 - 0.0281*log(dbh_cm)^3))
(AGB <- sum(BDplot.BD1$AGB)/1000) #Ton/ha
## [1] 469.1828
The ratio developed by Mokany et al. (2006). Critical analysis of root : shoot ratios in terrestrial biomes. Global Change Biology, 12: 84-96 offers specific ratios based on forest type and climate zone and are applicable when the aboveground biomass estimate (shoot) is reported at the stand level
(BGB = 0.235 * AGB)
## [1] 110.2579
Total biomass of 1 hecta
(TB <- AGB + BGB)
## [1] 579.4407
The quantity of carbon can then be estimated by converting biomass to carbon using the IPCC default carbon fraction of 0.47. After that, the carbon stock was multiplied by 44/12 (the carbon atom ratio in the molecular weight of CO2) to get the CO2 equivalent value.
#Carbon stock (ton/ha)
(TC <- TB * 0.47)
## [1] 272.3371
# CO2 (ton)
(CO2 <- TC * 44/22)
## [1] 544.6743
1 carbon credit = 1 ton of CO2 equivalent (tCO2e) emissions.
In 2023, Vietnam successfully sold 10.3 million forest carbon credits for the first time through the World Bank at a price of $5/ton (Global average selling price ~$11.2)
1 hecta forest at Bidoup - Nui Ba national park ~ 2700 $
(CO2 * 5)
## [1] 2723.371
LITHCA : Lithocarpus eucalyptifolius
CASTEC: Castanopsis echinocarpa
SYMPBA: Symplocos banaensis